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Finite-Difference  Time-Domain  Calculations 
Based  on  Recursive  Convolution  Approach 
For  Propagation  of  Electromagnetic  Waves 
in  Nonlinear  Dispersive  Media 

S.  J.  Yakura  and  J.  T.  McGillivray 
Air  Force  Research  Laboratory,  AFRL/DEPE 
Kirtland  AFB,  NM  87117-5776 
(30  October  1997) 

Abstract 

The  piecewise  linear  recursive  convolution  (PLRC)  approach  has  been  shown  to  provide 
much  improved  accuracy  over  the  usual  discrete  recursive  convolution  approach  while 
retaining  the  efficient  use  of  computer  memory  storage  and fast  computational  speed  for 
finite-difference  time-domain  electromagnetic  propagation  calculations  for  linear 
dispersive  materials.  In  this  paper,  an  idea  behind  the  implementation  of  the  PLRC 
approach  is  extended  to  handle  nonlinear  dispersive  media,  specifically  for  the 
convolution  integral  that  depends  on  the  product  of  the  electric  field  squared  and  the 
third-order  electric  susceptibility  function.  Compared  to  linear  dispersive  material, 
where  one  has  a  simple  linear  relationship  for  the  next  time  step  electric  field  as  a 
function  of  the  previous  time  step  electric  field,  the  nonlinear  dispersive  material  case 
has  a  cubic  equation  for  the  next  time  step  electric  field  as  a  function  of  the  previous  time 
step  electric  field.  Consequently,  the  cubic  equation  must  be  solved  at  successive  times  to 
advance  the  electric  field  in  the  next  time  step. 


I.  INTRODUCTION 

There  has  been  considerable  interest  in  understanding  the  transient  behavior  of  an 
ultrafast  laser  pulse  that  interacts  with  nonlinear  dispersive  materials.  In  the  last  several 
years  many  experimentalists  have  made  use  of  the  newly  developed  Kerr-lense  mode- 
locked  Titanium-Sapphire  lasers  to  perform  well-controlled  experiments  so  that  they  can 
obtain  accurate  transient  behavior  measurements  of  ultrafast  laser  pulses  in  simple 
molecular  liquids  and  solids  that  are  known  to  exhibit  nonlinear  optical  effects  [1].  To 
better  understand  the  details  of  the  nonlinear  processes  that  are  observed  in  experiments, 
numerical  simulation  has  been  used  extensively  in  reproducing  the  observed  nonlinear 
effects.  Until  recently  most  computer  simulation  has  been  performed  by  solving 
approximated  Maxwell’s  equations,  known  as  the  generalized  nonlinear  Schrodinger 
(GNLS)  equation  [2],  for  the  envelope  of  the  propagating  oscillating  wave  packet  that 
provides  information  on  the  time  evolution  of  the  overall  shape  of  the  optical  pulse.  Since 
computer  simulation  based  on  the  GNLS  equations  could  not  simply  provide  any 
information  about  oscillating  waves  contained  within  the  envelope  of  the  optical  pulse, 
there  is  a  renewed  interest  in  obtaining  the  greater  details  of  propagating  optical  pulse 
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phenomena  by  canying  out  computer  simulation  of  Maxwell’s  equations  directly. 

In  the  last  few  years,  computer  technology  has  advanced  to  the  point  where  arithmetic 
processing  chips  can  operate  up  to  hundreds  of  hertz  and  dynamic  random  memory  chips 
can  hold  in  excess  of  multi-giga  bytes  of  memory.  Using  the  present  day  computers,  we 
can  consider  solving  Maxwell’s  equations  directly  without  having  to  rely  on 
approximated  Maxwell’s  equations.  Of  recently  investigated  numerical  techniques  which 
showed  the  promising  future  is  the  well-known  finite-difference  time-domain  (FDTD) 
method  [3]  that  is  based  on  a  simple  differencing  scheme  in  both  time  and  space  to 
calculate  transient  behavior  of  electromagnetic  (EM)  field  quantities.  Because  of  the 
usefulness  of  the  FDTD  method  for  many  optical  applications,  recent  researchers  have 
focused  their  attention  into  numerical  handling  of  linear  and  nonlinear  polarization  terms 
which  appear  in  one  of  Maxwell’s  equations  as  convolution  integrals  so  that  they  can 
simulate  linear  and  nonlinear  dispersive  effects  more  effectively  [4-9]. 

Depending  on  the  form  of  the  integrand  appearing  in  the  convolution  integral,  the 
dispersive  effect  can  be  classified  as  linear  or  nonlinear.  For  linear  dispersive  materials, 
the  relationship  between  displacement  field  vector  U.  and  electric  field  vector  E  is  usually 
expressed  in  the  following  form 

D(t)  =  E(l)  +  ^, 

P  0 

where  €<,  is  the  electric  permittivity  of  free  space,  €„  is  the  permittivity  at  infinite 
frequency,  and  is  the  pth  term  first-order  electric  susceptibility  function  that 

depends  on  time  difference  (t-r). 

For  materials  that  show  both  the  linear  and  nonlinear  polarization  properties,  specifically 
through  the  first-order  (linear)  and  third-order  (nonlinear)  electric  susceptibility  functions, 
and  X^^^ p(t,x,ti,t 2),  respectively,  the  relationship  between  D  and  E  can  be 
expressed  in  the  following  form 

D(t)  = 

P  0 

+  JJ  )E.(f  )E.(f  ~f2 )x  p  (1'2) 

p  0  0  0 

where  X^^^ p(t,x,tj,t^  is  the  pth  term  four-time  dependent  third-order  susceptibility 
function  which  contributes  to  the  nonlinear  behavior  of  the  material. 

Reducing  X^^^p(t,z,tj,t2)  to  the  one-time  dependent  susceptibility  function,  X^^^p(trt2),  by 
use  of  the  following  Bom-Oppenheimer  approximation  [10], 

=  -t2)f  (1-3) 

where  aop  is  a  constant  and  b(t)  is  the  Dirac  delta  function,  we  can  reduce  Eq.(l  .2)  to  a 
more  amenable  expression  that  consists  of  the  sum  of  linear  and  nonlinear  convolution 
integrals  of  the  form 


2 


(1.4) 


t 

P  0 


*  S.  E(l)  2; 

P 


\o 


) 


Based  on  the  above  expression,  this  paper  provides  the  general  formulation  of  the  FDTD 
method,  which  we  call  the  piecewise  continuous  recursive  convolution  (PCRC)  method, 
to  evaluate  the  linear  and  nonlinear  single  convolution  integrals.  We  investigate 
specifically  the  case  in  which  both  the  first-order  and  third-order  electric  susceptibility 
functions  are  expressed  in  the  following  exponential  forms 


X  (t)  =  Real  {apexp[^  p  t]}U(t)  (1.5) 

x[^^(t)  =  Real  {al‘expH’';t]}U(t)  (1.6) 

where  U(t)  is  the  unit  step  function  and  a^p,  a"‘p,  y^p  and  y"  p  are  complex  constants. 

By  making  proper  choices  of  complex  constants  and  carrying  out  the  Fourier  Transform, 
we  can  readily  obtain  the  familiar  Debye  and  Lorentz  forms  of  the  complex  permittivity 
from  the  above  exponential  forms. 


II.  GOVERNING  EQUATIONS  AND  GENERAL  FORMULATION 

Starting  with  Maxwell’s  equations,  we  can  express  the  following  differential  equations 
for  spatial  and  time  dependent  EM  field  quantities  inside  the  dispersive  material. 


YxH(t)= 


Y^E(t)  =  - 


dt 

d(\iH(t)) 

dt 


E(t)  +  e„  Yp  (t) 


+  E(t)Y^  pyt)+ail^j[E(x)f5(t-x)dx 

P  V  0 

0 

pyt)J\[E(x)rx<;>(t-x)dx 


(2.1) 

(2.2) 

(2.3) 

(2.4) 

(2.5) 
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where  H  is  the  magnetic  field  vector,  D  is  the  displacement  field  vector,  E  is  the  electric 
field  vector,  a  is  the  electrical  conductivity,  p  is  the  magnetic  permeability,  Gq  is  the  free 
space  electric  permittivity,  e«>is  the  permittivity  at  infinite  frequency,  and  are 
the  pth  term  first-order  and  third-order  electrical  susceptibilities,  and  P  p  and  P^'p  are 
related  to  the  pth  term  linear  and  nonlinear  polarization  field  vectors. 

Using  FDTD  the  above  equations  can  be  solved  numerically  at  each  time  step,  provided 
P’p(t)  and  P"‘p(t)  are  handled  properly.  Thus,  the  whole  problem  rests  upon  proper 
numerical  evaluation  of  P’p(t)  and  P"^p(t)  at  each  time  step.  For  that  reason,  the  rest  of  this 
section  is  devoted  to  the  discussion  of  the  numerical  formulation  used  to  solve  P'p(t)  and 

In  order  to  achieve  better  accuracy  in  evaluating  the  convolution  integral.  Eft)  is 
considered  to  be  a  piecewise  continues  function  over  the  entire  integration  limits  and  Eft) 
changes  linearly  with  respect  to  time  over  a  given  discrete  time  interval  [wAf, 
where  m  represents  the  previous  mth  time  step  of  the  current  nth  time  step  in  the  FDTD 
calculation  [1 1].  Referring  to  Fig.l,  in  terms  of  the  electric  field  values,  and  Er^\ 
which  are  evaluated  at  discrete  time  steps  t=mAt  and  t=fm+l)At,  respectively,  we  can 
express  Eft)  in  the  following  form. 


Eft)  =  E""' 


fE"‘^‘-E”') 

At 


ft-mAt); 


mAt<t  <  f  m  +  l)At  <  nAt 


(2.6) 


When  Eq.(2.6)  is  substituted  into  Eq.(2.4),  we  can  obtain  the  following  expressions  for 
discrete  values  of  fP^p}"  at  time  increment  n  (see  Appendix  for  manipulations). 

(K)‘=Y.{E‘'('ii‘,,r-[e"‘-E’](V,jr}  (2.7) 

m=0 


where 

(n-m)At 

\  X'p'YtM  (2.8) 

j  (n-m)td 

J  [(r>-m)LI-x]X<>(t)ch  (2.9) 

Similarly,  substituting  Eq.(2.6)  into  Eq.(2.5),  we  can  obtain  the  following  expressions  for 
discrete  values  of  (P”y"  at  time  increment  n  (see  also  Appendix  for  manipulations). 

(pir  =i^{ 

m=0 

-2  -f:'”)  (xi/p^, )'*’'"  --£^'”)^(m/^2)"’'”}  (2.10) 

where 
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(2.11) 


J 


(n-l-mJAi 


(  n-m 


(2.12) 


1 


(n-m  )At 

J  [(n-m)^t-x]^  X^^^(x)dx 

(n-l-m  )At  ^ 


(2.13) 


When  Eq.(1.5)  is  substituted  in  Eqs.(2.8)  and  (2.9),  and  Eq.(1.6)  in  Eqs.(2.1 1),  (2.12) 
and  (2.13),  we  can  show  after  some  manipulations  the  following  recursive  relationships  to 
exist  between  next  time  step  («+i)Ar  and  current  time  step  nAt. 


(yvlo)""'''”  =exp(-ylAt)(xvU'” 

(2.14) 

=exp(-ylAt)(yvl,)"'” 

(2.15) 

=exp(-y;^At)(^,;^X"' 

(2.16) 

rvF”pV”^^'”=ex;,r-yp"^Ar;rvK;V”'” 

(2.17) 

^exp(-y;^At)(x^l,r 

(2.18) 

We  are  able  to  obtain  the  above  recursive  relationships  only  because  the  susceptibility 
functions  are  expressed  in  the  exponential  form. 

When  Eqs.(2.14)  through  (2.18)  are  used  in  Eqs.(2.7)  and  (2.10),  respectively,  for  the 
next  discrete  time  step  at  (n+l)At,  we  can  obtain  the  following  recursive  relationships  for 
and  in  terms  of  and 


(p[r‘ 


= -  E”)(^ 
+  exp(-y'^At)(p^p)" 


(2.19) 


=rr"'/('v"pV'’  -  2E’’*‘(E”*‘ -E'‘)(y!j) 


+ (E”^‘  -  EV'  (yv : V"  +  ^^P(-fo  A 0  (Pi') 


nt  \n-k-l,n 
P.^ 

-,ni  \n 


(2.20) 


In  the  above  expressions,  (¥”^0)”"^’"  (V” and 

are  evaluated  at  each  succesuve  time  step  using  the  recursive  relationships  [see  Eqs.(2.14) 
through  p.  18)]  starting  with  initial  values  of  (Vp_o)^'"  (v|/p,/)^’"  (v”p,o)^’”  (ij/"  and 
(V"  p,2)  ’  >  which  are  calculated  explicitly  at  the  beginning  of  computer  simulation  for  the 


5 


selected  value  of  At  [see  Eqs.(2.8),  (2.9),  (2.1 1),  (2.12)  and  (2.13)  and  set  n=l  and  m=0]. 
To  demonstrate  how  the  above  terms  are  used  in  FDTD  calculations,  we  considered  the 
one  dimensional  case  where  the  propagating  wave  vector  k  lies  along  a  major  axis  in 
Cartesian  coordinates.  The  analysis  can  be  extended  easily  in  the  three  dimensional  case 
where  the  EM  wave  is  propagating  in  any  arbitrary  direction  based  on  the  following 
sample  formulation  as  described  below.  For  our  one  dimensional  case,  we  arbitrarily 
picked  the  k  vector  to  lie  in  the  x  direction,  Dy  and  Ey  field  vectors  in  the  y  direction  and 
field  vector  in  the  z  direction.  When  Eq.(2. 1)  is  differenced  in  both  time  and  space 
using  the  usual  Yee  algorithm,  we  have 

At  Ax  2  {2.21) 

where  n  and  i  indices  are  used  to  denote  discrete  nth  time  step  nAt  and  ith  spatial  location 
iAx,  respectively. 

Using  Eqs.(2,3),  (2.19)  and  (2.20),  the  left-hand  side  of  Eq.(2.21)  can  be  expressed  in 
terms  of  and  (ETjj  as  follow 

+  e.  Safi;*' -[(E";')-(E;),](^ur''’ 

P 

^[exp(-ylM)-I][(p',r],  } 

[<ET‘),r  -  [(e;),]’j  aj'p'  +  re;'’ ),r(v:‘,r‘- 

p 

-  2  [(e;^^  r  [(E-;' -  (e;  m  t,  r'" 

+[(E”;‘hexp(-^”;At) -(E;)j[(Pi^r], }  (2.22) 

When  Eq.(2.22)  is  substituted  in  Eq.(2.21),  it  results  in  the  following  cubic  equation  in 
which  we  need  to  solve  for  {Ey*‘)i . 

2  o. =  0  (2-23) 

k=0 

where 

<HT’‘ +  <e;),^-  (E;),s,e. 

Ax  2 

P 

^  i(E‘),r  ]  «!.'>  +(■£;;, /'rr// 7,;  (2.24) 
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CT 


a,=- 


+€„g„+£„2;  ((^lor''’  -  (wiir''"  + 

p 

+  exp(-y;‘AOf(Pl'rJ,  ) 


<-2 'L( 2(E:),<w“,r''' -  2(E;M’v;‘,r''} 

p 

=s.  'Ef  ai?+  rv".r" -2(<if;‘,r‘'  +('vf.2r‘') 


(2.25) 

(2.26) 

(2.27) 


Note  that  because  we  defined  the  electric  susceptibility  functions  to  be  real  parts  of 
complex  exponential  functions  as  shown  in  Eqs.(1.5)  and  (1.6),  it  is  important  that  we 
choose  only  the  real  parts  when  evaluating  the  above  constants  ao,  aj,  a2,  and  a.^. 

The  above  equation  can  be  solved  for  (Ey^^)i  using  any  one  of  the  root  finding 
numerical  techniques. 

At  each  successive  time  step,  only  Eqs.(2.19),  (2,20)  and  (2.23)  have  to  be  solved  for 
updated  /■(P^p)”'^^7j,  and  values  in  order  to  handle  the  electric  field 

response  of  nonlinear  dispersive  materials. 

For  the  purely  linear  dispersive  case,  a2  and  a3  as  well  as  some  terms  appearing  in  a^  and 
ai_  turn  out  to  be  zero.  In  this  case  we  can  solve  for  (Ey^‘)i  directly  without  having  to  rely 
on  the  numerical  root  finding  technique  as  seen  in  many  previously  published  papers 
[12-18]  that  discuss  computational  schemes  for  linear  dispersive  materials. 


III.  NUMERICAL  ANALYSIS  -  A  CASE  STUDY  FOR  NONLINEAR  SOLITON 
FORMATION 

To  demonstrate  the  validity  of  the  PCRC  method,  we  investigated  the  case  of  an  optical 
pulse  where  it  propagated  in  the  x-direction  in  free  space  and  incident  on  an  infinite  half 
space  dispersive  medium  that  is  characterized  by  zero  electrical  conductivity  and  the 
following  single  time  dependent  first-order  (linear)  and  third-order  (nonlinear) 
susceptibility  functions,  and  p(t)  [19]. 

For  linear  dispersion  contribution: 

y^(i)(t)J.^Al-exp(^t)sin(v„t)  (3.1) 

^  o 

For  nonlinear  dispersion  contribution  arising  from  purely  Raman  scattering  and  no 
virtual  electronics  transition  effect  [20]; 

=Xo^ [(U  *A)hi^l]exp(-t/x2)sin(t  jxj)  and  a[^l  =  0  (3.2) 

where  ©/s  e  Vo^Vfioo^-S^),  ©^  is  the  resonant  frequency,  is  relative 
permittivity  at  DC,  6  is  the  first-order  susceptibility  damping  constant,  is  the 
nonlinear  coefficient,  7/t;  is  the  optical  phonon  frequency,  and  X2  is  the  optical  phonon 
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lifetime. 

Comparing  Eqs.(3.1)  and  (3.2)  with  the  previously  defined  susceptibility  functions  as 
seen  in  Eqs.(1.5)  and  (1.6),  we  can  relate  the  above  coefficients  to  the  previously  defined 
susceptibility  coefficients  as  follows. 


I 

ap 


m  \/v  „  and 


Yp 


-("5  +/vj, 


ap  [(A 


(3.3) 

(3.4) 


where  i  is  the  imaginary  number  defined  as  V-1 . 

The  incident  optical  pulse  is  assmned  to  propagate  at  the  sinusoidal-carrier  electric  field 
frequency,  of  8.61x10*'*  rad/sec  with  the  unit  amplitude  (1.0  volts/meter),  enveloped 
inside  the  hyperbolic  secant  function  represented  by  characteristic  time  constant  T^.  is 
the  parameter  used  to  determine  the  width  of  the  overall  shape  of  the  incident  optical 
pulse.  Thus,  we  used  the  following  expression  to  launch  the  incident  optical  pulse  in  free 
space  for  our  FDTD  simulation. 


Incident  Optical  Pulse  (t)  =  cos  [<o^  (t-  tj^,^  ) ]  sech\ 


( I  ~  I  delay  ) 

T 


(3.5) 


where  is  the  delay  time  for  the  incident  optical  pulse  to  reach  the  peak  value  at  the 
place  where  the  optical  pulse  is  laimched.  Two  different  values  of  the  characteristic  time 
constant  are  used  to  investigate  the  effect  of  the  characteristic  time  constant  on  soliton 
formation.  We  used  3.5  femtoseconds  (fs)  and  7.0  fs  for  the  characteristic  time  constants. 
These  characteristic  time  constants  resulted  in  enveloping  about  6  cycles  of  the  electric 
field  oscillation  for  the  3.5  fs  pulse  and  12  cycles  for  the  7.0  fs  pulse. 

We  selected  the  total  simulation  cells  to  be  50,000,  ranging  from  x=- 10,000  to  x=40,000 
with  the  fi'ee  space-dispersive  material  interface  located  at  x=0.  The  optical  pulse  was 
launched  in  firee  space  at  x=- 10,000,  traveling  in  the  positive  x  direction.  The  LIAO 
absorbing  boundary  condition  [21]  was  used  at  the  two  end  points  of  the  computational 
space. 

For  the  selection  of  basic  FDTD  parameters  we  chose  the  following  values. 

Uniform  cell  size  (or  Ax)  =  5.0  nanometers  (  =  2'Kc/GiJA'h%  =  X/438), 

“18 

Time  step  increment  (or  At)  =  8.33  attoseconds  (10'  )  ( =  hx/2c\ 

Total  number  of  uniform  cells  =  50,000, 

Total  number  of  time  steps  =  2x10^, 

where  c  is  the  speed  of  light.  For  Ax  of  5  nanometers,  we  can  estimate  the  free  space 
numerical  phase  velocity  error  to  be  aroimd  5x10'^  [22]. 

To  observe  the  soliton  formation  within  the  total  propagation  distance  of  250 
micrometers  (=5.0  nanometers/cell  x  50,000  cells),  we  had  to  enhance  the  nonlinear 
behavior  by  scaling  the  coefficients  foimd  in  the  first  and  third  order  susceptibility 
functions  (see  Eqs.(3.1)  and  (3.2)).  To  show  consistencies  of  our  FDTD  results  with 
previously  published  results,  we  used  the  values  similar  to  the  ones  found  in  the  papers 
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referenced  in  [4,5].  Shown  below  are  the  values  that  we  used  for  our  FDTD  calculations. 

For  linear  dispersive  material  properties: 

=  5.25;  €„  =  2.25; 

©o  ==  8x10^'*  see'*;  6  =  4.0x10^  sec'*. 

For  nonlinear  dispersive  material  properties: 

=  30.0  (volts/meter)^; 

T;  =  12.2  fs;  T2  =  32.0  fs. 

To  see  the  difference  in  linear  and  nonlinear  responses,  we  first  calculated  the  linear 
dispersive  case  by  setting  the  nonlinear  coefficient,  to  zero.  For  the  nonlinear  case, 
we  used  the  nonlinear  material  property  values  as  shown  above. 

To  calculate  the  updated  electric  field  value  at  each  time  step,  we  used  the  simple 
Newton  iteration  method  by  making  use  of  the  previous  time  step  electric  field  value  as 
the  initial  guess  to  solve  the  cubic  equation  (see  Eq.(2.23)).  For  every  simulation  we 
performed,  the  convergence  criterion  of  lO"^  was  reached  after  at  most  two  iterations. 

Shown  in  Fig.2  and  Fig.3  are  spatial  plots  obtained  for  linear  and  nonlinear  dispersive 
calculations,  respectively,  at  time  steps  of  10,000, 40,000, 70,000, 100,000  and  130,000 
for  the  incident  optical  pulse  width  of  3.5  fs.  Similarly,  Fig.4  and  Fig.5  show  spatial  plots 
obtained  for  the  same  input  values  with  the  exception  of  7.0  fs,  instead  of  3.5  fs,  for  the 
incident  optical  pulse  width.  As  seen  in  Fig.2  and  Fig.3  (likewise  Fig.4  and  Fig.5)  the 
solitary  wave  started  forming  first  with  the  appearance  of  the  small  spike-like  shape 
inside  the  linearly  dispersive  part  of  the  pulse  as  a  result  of  the  nonlinear  self-focusing 
effect.  As  the  pulse  propagated  deeper  in  the  dispersive  medium,  the  spike-like  shape 
transformed  gradually  to  the  shape  that  resembles  the  solitary  wave  packet  and  became 
isolated  from  the  main  linear  dispersive  part  of  the  pulse  due  to  the  slower  phase  velocity. 
Once  the  solitary  wave  packet  became  completely  isolated  from  the  linear  dispersive  part 
of  the  pulse,  the  solitary  wave  packet  propagated  at  constant  amplitude  while  maintaining 
the  general  structure.  On  the  other  hand,  the  linear  dispersive  part  of  the  pulse  decreased 
in  its  amplitude  and  became  much  broader  in  its  shape  as  it  propagated  deeper  into  the 
dispersive  medium  because  of  the  linear  dispersive  effect. 

Shown  Fig.6  is  the  comparison  of  two  solitary  wave  packets  that  are  formed  from  the 
two  different  incident  optical  pulse  widths.  It  shows  the  overlay  views  taken  around 
solitary  wave  packets  of  the  two  spatial  plots  which  are  found  in  Fig.3  and  Fig.5 
(specifically  at  the  time  step  of  100,000).  We  can  see  that  these  two  solitary  wave  packets 
have  approximately  the  same  size  envelopes.  We  estimated  the  solitary  wave  packet  sizes 
to  be  arormd  1.54  micrometers  based  on  approximately  two  wave  cycles  of  oscillation  at 
a  wavelength  of  0.77  micrometers  contained  inside  the  envelopes  of  these  solitary  wave 
packets.  We  can  think  of  that  the  formation  of  the  same  size  solitary  wave  packet  is 
analogous  to  calculating  the  allowable  bound  states  of  the  nonlinear  Schrodinger  wave 
equation  inside  a  potential  well  [23].  Because  of  material  property  values  used  for  our 
simulation  calculations,  we  ended  up  getting  the  unique  solitary  wave  packet  size  as  we 
saw  in  our  simulation  results. 

Also  shown  in  Fig.6  is  the  difference  we  obtained  in  the  amplitude  of  the  two  solitary 
wave  packets  from  two  different  incident  optical  pulse  widths.  The  wider  incident  pulse 
resulted  in  about  1 .2  times  that  of  the  narrow  incident  pulse. 
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For  the  nonlinear  case,  we  saw  the  formation  of  a  small  secondary  high  frequency  pulse 
that  moved  ahead  of  both  the  linear  dispersive  part  of  the  pulse  and  the  nonlinear  solitary 
wave  packet.  The  same  secondary  high  frequency  pulse  was  also  obtained  and  reported 
previously  for  computational  modeling  of  femtosecond  optical  solitons  using  another 
FDTD  approach,  called  the  auxiliary  differential  equation  approach,  to  handle  the 
nonlinear  dispersive  term. 

To  perform  these  calculations  we  used  a  SPARC  20  workstation  equipped  with  a  75 
MHz  processor  and  512  Megabyte  dynamic  random  access  memory  chips.  On  average  it 
took  around  24  CPU  hours  to  complete  the  job  with  no  optimization.  There  was  little 
difference  in  the  total  computational  time  for  the  linear  dispersive  case  as  compared  to  the 
nonlinear  dispersive  case.  Using  present  day  computers  equipped  with  much  faster 
multiple  processors,  the  computational  time  can  be  reduced  considerably  more. 


IV.  CONCLUSIONS 

In  conclusion,  we  have  shown  in  our  sample  calculations  that  the  PCRC  approach  of  the 
FDTD  method  presented  here  is  fully  capable  of  predicting  the  formation  of  nonlinear 
solitary  waves  by  solving  Maxwell’s  equations  directly.  The  PCRC  approach  resulted  in  a 
much  simpler  algebraic  form  to  relate  the  displacement  field  vector  to  the  electric  field 
vector  than  the  auxiliary  differential  equation  approach  which  requires  the  additional 
coupled  nonlinear  ordinary  equations  to  be  solved  at  each  time  step.  Consequently,  the 
PCRC  approach  is  capable  of  calculating  at  much  faster  computational  speed.  Also, 
because  of  the  piecewise  linear  approximation  used  for  the  time  dependent  electric  field 
vector,  the  PCRC  approach  should  provide  the  accuracy  comparable  to  that  of  the 
auxiliary  differential  equation  approach. 

Also,  we  gained  much  from  exponential  function  forms  of  the  linear  and  nonlinear 
susceptibility  functions  which  allowed  us  to  implement  the  recursive  feature  in  our 
algorithm.  As  a  whole,  the  PCRC  approach  retained  all  the  advantages  of  the  usual 
discrete  recursive  convolution  approach,  such  as  fast  computational  speed  and  efficient 
use  of  the  computer  memory,  however,  with  much  improved  accuracy. 

The  nonlinear  dispersive  formulation  resulted  in  having  to  solve  the  cubic  equation  for 
the  successive  time  step  electric  field  values  as  compared  to  the  linear  equation  for  the 
simple  linear  dispersive  case.  For  sample  calculations  we  have  looked  into  here,  the 
simple  Newton  iterative  method  provided  sufficiently  fast  convergence  for  finding  the 
root  of  the  cubic  equation  by  using  the  previous  time  step  electric  field  value  as  the  initial 
guess. 

The  present  PCRC  approach  is  robust  and  applicable  for  applications  in  two  and  three 
dimensional  problems.  The  one  dimensional  code  can  be  extended  easily  into  two  and 
three  dimensional  codes  with  little  effort. 
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Fig.  1 .  Illustration  of  Piecewise  Linear  Approximation  for  the 
Electric  Field  as  a  Function  of  Discrete  Time  Steps 
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Fig,  2.  Perspective  plots  of  the  electric  field  intensity  versus  spatial  location  from  linear 
dispersive  calculations  taken  at  five  successive  times  of  Ti,  T2,  T3,  T4  ant  T5  in 
order  to  show  the  space-time  evolution  of  the  optical  pulse  with  the  initial  pulse 
width  of  3.5  fs. 

[  Ti=  10,000At  (=0.083  ps);  T2=  40,000At  (=0.333  ps); 

T3=  70,000At  (=0.583  ps);  14=100, OOOAt  (=0.833  ps); 

T5=130,000At  (=1.083  ps)  with  At=8.33xl0'^*  second  ] 
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Electric  Field  (v/in) 


Fig.  3.  Perspective  plots  of  the  electric  field  intensity  versus  spatial  location  from 

nonlinear  dispersive  calculations  taken  at  five  successive  times  of  Ti,  T2,  T3,  T4 
ant  T5  in  order  to  show  the  space-time  evolution  of  the  optical  pulse  with  the 
initial  pulse  width  of  3.5  fs. 

[  Ti=  10,000At  (=0.083  ps);  T2=  40,000At  (=0.333  ps); 

T3=  70,000At  (=0.583  ps);  14=100, OOOAt  (=0.833  ps); 

T5=130,000At  (=1.083  ps)  with  At=8.33xl0’^*  second  ] 
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Electric  Field  (v/m) 


Fig.  4.  Perspective  plots  of  the  electric  field  intensity  versus  spatial  location  from  linear 
dispersive  calculations  taken  at  five  successive  times  of  Ti,  T2,  T3,  T4  ant  T5  in 
order  to  show  the  space-time  evolution  of  the  optical  pulse  with  the  initial  pulse 
width  of  7.0  fs. 

[  Ti=  10,000At  (=0.083  ps);  T2=  40,000At  (=0.333  ps); 

T3=  70,000At  (=0.583  ps);  T4=100,000At  (=0.833  ps); 

Ts  =1 30, OOOAt  (=1.083  ps)  with  At=8.33xl0'**  second  ] 
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Fig.  5.  Perspective  plots  of  the  electric  field  intensity  versus  spatial  location  from 

nonlinear  dispersive  calculations  taken  at  five  successive  times  of  Tj,  T2,  T3,  T4 
ant  Ts  in  order  to  show  the  space-time  evolution  of  the  optical  pulse  with  the 
initial  width  of  7.0  fs. 

[  T,=  10,000At  (=0.083  ps);  T2=  4O,000At  (=0.333  ps); 

T3=  70,000At  (=0.583  ps);  T4=100,000At  (=0.833  ps); 

Ts  =130, OOOAt  (=1.083  ps)  with  At=8.33xl0'^*  second  ] 
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Electric  Field  (v/m) 


Fig.  6.  Comparison  of  solitary  wave  packets  obtained  in  the  nonlinear  dispersive  medium 
from  launching  two  different  width  pulses  in  free  space.  The  solid  line  is  from  the 
optical  pulse  width  of  7.0  fs  and  the  dashed  line  is  from  the  optical  pulse  width  of 
3.5  fs.  These  are  taken  at  the  time  step  of  100,000. 
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APPENDIX 


For  the  linear  part,  when  we  substitute  Eq.(2.6)  in  Eq.(2.4)  and  evaluate  t  at  nAt  for  the 
argument  of  the  susceptibility  function  where  n  is  the  nth  time  step,  we  have 


n-1 

m=0 


m+1 


-E’")^^[x-mAt]  }X  l'^(nAt-x)dx 


n-l  (m+l)ls.t 

=  Y.(E'”*'  I  Xl'^(nAt-x)dx 

m=0 


-  -E'")—  j  [x-mAt]X['^(nAt-x)dx} 

mAl 


(A.1) 


Similarly  for  the  nonlinear  part,  when  we  substitute  Eq.(2.6)  in  Eq.(2.5)  and  evaluate  t  at 
nAt  for  the  argument  of  the  susceptibility  function  where  n  is  the  nth  time  step,  we  have 


M-/  ( J 

( =  S  W (E”'^‘f  -  2  E”'*’  (E'”^’  -E”')  — 

■n=0  miJ  ^  * 


+  ^  [x -mAt]' }x['^(nAt-x )dx 


n-1 


(  m+1  )Ai 

=  I  xi”<'<M-r)dx 

m=0  wA/ 


-2E”‘*'(E”'^‘ -E’”)—  \  [x-mAt]'i<J^(nAt-x)dx 


(  m+l  )At 


At 


mtJ 


+(e”'^'-e'”  y  —E 


(m+})At 


(At )  mAt 


J  [x -mAt]^  Xp^(nAt-x  )dx  } 


(A.2) 


Using  the  change  of  variable  x  ’=(nAt-x),  we  can  readily  show  the  existence  of  the 
following  relationships  for  the  above  integrals 

(m+1)  At  (n-m)At 

J  [x -mAt]^  f  (nAt -X  )dx  =  j  [(n-m)At-x'  f(x'  )dx'  (A.3) 

mAt  (n-l-m)At 


where  k  takes  the  values  of  0, 1,  and  2  and  f(t)  represents,  respectively,  the  time- 
dependent  first  and  third  order  susceptibility  functions,  and 

Since  t  ’  appearing  in  the  right-hand  side  of  Eq.(A.3)  is  the  integration  variable,  we  can 
simply  replace  it  by  x.  When  Eq.(A.3)  is  substituted  back  in  Eqs.(A.l)  and  (A.2),  we  can 
obtain  the  expressions  as  shown  in  Eqs.(2.7)  and  (2.10). 
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